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A large class of quantum and statistical field theoretical models, encompassing relevant condensed 
matter and non-abelian gauge systems, are defined in terms of complex actions. As the ordinary 
Monte-Carlo methods are useless in dealing with these models, alternative computational strategies 
have been proposed along the years. The Langevin technique, in particular, is known to be frequently 
plagued with difficulties such as strong numerical instabilities or subtle ergodic behavior. Regarding 
the chirally decomposed version of the sine-Gordon model as a prototypical case for the failure of the 
Langevin approach, we devise a truncation prescription in the stochastic differential equations which 
yields numerical stability and is assumed not to spoil the Berezinskii-Kosterlitz-Thouless transition. 
This conjecture is supported by a finite size scaling analysis, whereby a massive phase ending at a 
line of critical points is clearly observed for the truncated stochastic model. 
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I. INTRODUCTION 

Several field theory models devoted to concrete con- 
densed matter and high energy applications are given 
in terms of complex- valued bosonic ("c- number") ac- 
tions. In the condensed matter context, for instance, 
they are the focus of much attention in low-dimensional 
systems related to Luttinger liquids and quantum Hall 
phenomena T, "5, "s", . On the high-energy front, on the 
other hand, a fundamental open issue refers to the phase 
diagram of SU(3)-QCD with non-zero fermion density, 
where the gauge action notoriously acquires an imaginary 
part, due to the non-vanishing quark chemical potential 

iiiiiimiia. 

Numerical non-perturbative efforts to study all of the 
above problems are considerably restricted by the fact 
that the Monte-Carlo method cannot be straightfor- 
wardly applied, since the weights in the path-integration 
measures are not positive-definite quantities anymore. 
To overcome this difficulty, a number of alternative nu- 
merical methods have been suggested, from time to time, 
such as complex Langevin simulations Il2l Il3l |. analyti- 
cal expansions [l^ . kernel mappings |15|. microcanon- 
ical averagin g pr odecures 16], path-integral factoriza- 
tion schemes [13, etc. To date, there are not yet well- 
established results on how general and efficient these 
tools are. Among them, the Langevin approach - the 
central subject of discussion in this paper - has as at- 
tractive features its great simplicity and direct numerical 
implementation. However, as it has been learned from 
purely empirical investigations, the stochastic evolution 
provided by the Langevin equations is affected, in many 
important cases, by bad numerical convergence and the 
existence of blow-up trajectories 0,(T8'|. Also, some spe- 
cial care may be needed when writing the Langevin equa- 
tions for actions which have cuts in the complex plane, 
in order to assure good ergodic properties of the phase 
space flow 0,13. 

It is not hard to understand in an essential way the 
origin of eventual instabilities in complex Langevin equa- 



tions. Since the field variables are promoted to complex 
numbers in the stochastic differential equations, the con- 
figurational phase space is enlarged, and unstable man- 
ifolds appear along the imaginary direction, containing 
the classical vacuum solutions. The time evolution of 
a phase space point (that is, a field configuration) will 
depart, in brownian-like motion, from the determinis- 
tic trajectory that would be followed in the absence of 
the stochastic noise. Thus, as soon as the phase space 
point happens to visit a region of "intense stream" in 
the neighborhood of an unstable manifold, the numerical 
convergence of the Langevin equation is under risk. 

While in polynomial models like the A^^ with complex 
coupling constants the Langevin simulation is mostly suc- 
cessful in the cases of physical interest, difficulties arise 
in other well-known models where the self-interaction po- 
tential is non-polynomial and invariant under compact 
Lie groups. The Langevin equations will contain, due to 
the extension of fields to the complex plane, "explosive" 
hyperbolic functions, leading, in general, to unstable nu- 
merical simulations 18J. That is precisely the situation 
found in many of the bosonized condensed matter sys- 
tems or in the finite Baryon density QCD realm. 

In order to investigate the problem of numerical 
Langevin stability, we have regarded the exactly solv- 
able euclidean sine-Gordon (SG) model as an ideal test- 
ing ground prototype, aiming, in this simpler context, 
to antecipate applications for a large class of systems. 
As briefly reviewed in Sec. II, the SG partition func- 
tion is alternatively written in terms of interacting left 
and right chiral bosons, governed by a complex- valued 
action. In Sec. Ill, we introduce the stochastic quanti- 
zation procedure for the chirally decomposed SG model, 
which is observed to produce numerical blow-up solu- 
tions. It turns out that the instabilities are not fixed by 
a reasonable time step resizing, or even with the help 
of more precise evolution algorithms (as a fourth-order 
Rungc-Kutta scheme). We propose, then, a pragmatical 
solution of the numerical convergence problem, through a 
straightforward "truncation prescription" applied to the 
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stochastic equations. The basic idea is to modify the 
Langevin equations without spoihng its thermodynami- 
cal phases. In fact, as a strong piece of evidence sup- 
porting our method, the Berezinskii-Kosterhtz-Thouless 
(BKT) phase transition pol I2H of the original SG model 
is verified by means of a finite size scaling study. Finally, 
in Sec. IV, we comment on our findings and discuss di- 
rections of further research. 



II. THE CHIRALLY DECOMPOSED SG MODEL 

As it is well-known, the usual euclidean SG correlation 
functions may be obtained, in path-integral language, 
from the zero temperature partition function p^ . 



Z = J D^cM-J d'xi^id^)^ _ |cos(/30)]} . (2. 



1) 



It is worth summarizing here some important facts on 
the above model |2^. In the g ^ limit, the SG spec- 
trum is massless for > Stt = /?^, defining a line of 
critical points, and massive for < (3c- f3 = f3c the 
SG model has a BKT phase transition. The same behav- 
ior holds for finite g, with /3c = Pc{g) (/3c grows with g). 
The SG model may be regarded as the bosonized version 
of the massless Thirring model perturbed by a Gross- 
Neveu interaction, defined in terms of a Dirac fermion 
field tp Its phases can be distinguished by the mass 
order parameter m = (i^ip), which is proportional to 
<I> = (cos(/3(/)/2)). For /3 — > /3c from below, the exact 
solution of the SG model yields ('0V') ~ exp[c/(/3 — /3c)] . 
For /3 > /3c the SG infrared fixed point is just the gaussian 
model and the order parameter vanishes in the thermody- 
namical limit. At /3 = /3c the order parameter has scaling 
dimension A(<I>) — —1/2, that is, for a finite size system 
with linear dimension L, it follows that ($) L~^/^. 

The chirally decomposed version of the SG model 
(CDSG) can be readly derived in a few elementary steps. 
Introducing the conjugate momentum tt = 7r(a;o,a:i), we 
may rewrite ()2.1() as 

Z = j D(I)Dt: exp{~ J (f + indocj) 
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+ ^(a<A)^-|cos(/3<^)]} . 
Define now the right and left chiral fields. 



(2.2) 
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C^f7I"(xo,0] 



(2.3) 



It is straightforward to invert H2.3|) . We find 

4> = 4>R + 4'L , 

TT = di(l)R - di(i)L ■ 



(2.4) 



As the jacobian of the mapping {4),i:) —>■ {4'Rt4'l) does 
not depend on the fields, the partition function (|2.2() may 
be expressed, by means of a direct substitution, as 



Z = j D(t>RD4>LGM-S) , 



(2.5) 



where 



+ di4>L{di - ido)^L - I cos(/3(0ij -I- </)l))] (2.6) 

is, thus, defined as the action of the CDSG model. It is 
clear that correlation functions of the SG model can be 
computed, in principle, from the above partition func- 
tion, through the substitution (f> — > <j)R + 4>l- Observe, 
however, that the CDSG action is complex, and the nu- 
merical chiral correlation functions cannot be obtained 
with the help of ordinary Monte-Carlo techniques. As 
an alternative numerical method, we will investigate the 
CDSG model, in the next section, from the point of view 
of Langevin equations, based on the formalism of stochas- 
tic quantization. 



III. STOCHASTIC SIMULATIONS 

The Parisi-Wu stochastic quantization method [25l |. 
states, under very general hypothesis, that euclidean field 
theories can be modelled through Langevin {— stochas- 
tic) differential equations. Actually, the CDSG model is 
equivalently given by the following set of Langevin equa- 
tions: 

OrCPR = -T— +VR , 
0<PR 

Ot<Pl = -7-7- +VL , 
OCPL 



(3.1) 



where r is the auxiliary "stochastic time" . In the asymp- 
totic T —> 00 limit, we expect that correlation functions 
involving the chiral fields, 4>r(t, xq, xi) and (/)l(t, xq, xi), 
taken at equal stochastic times, will converge to the 
CDSG ones. Both tjr and 77^ are gaussian random fields, 
satisfying to 

{iIr) ^ (vl) = , 
{VRiT,x)T]LiT',x')) = , 
{VR{r,x)7jRiT',x'))^2S{r-T')S\x-x') , 
{VLir,x)7jL{T\x'))=2S{T-T')S'{x~x') . (3.2) 

Above, X and x' denote points in the physical two- 
dimensional space. Substituting the action H2.6(l into 
lprT|l we get 

dr(l)R = 2di{di + idQ)(j)R - 5 sin(/30) + rjR , 
dr(l)L = 2di{di - idQ)(j)L - gsin(/30) -I- tjl , (3.3) 
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where 4> = 4'R~^4'l- Since these equations are complex, it 
is convenient, for simulation purposes, to split their real 



and imaginary components. Define 4>r 



dr' 



+ i> 



4^^ 



i4>'^ and 

The resulting equations are, therefore, 

5Cosh(/30(2))sin(/30(i)) + ryH , 

gsinh(/3(^(2))pQg(^^(i)) ^ 

292^^1) + 2didocl)f 
gcosh(/3(^(2)) gjj^(^0(i)) _,_ ^ 

2dl<i)f + 2dido(l>^l^ 
- gsinh(/3(/)(^))cos(/30(^^) . (3.4) 

It is a hopeless task to solve the above coupled differ- 
ential equations via direct numerical simulations. In gen- 
eral, for /3 close to [3c (but not necessarily within the crit- 
ical region), the numerical convergence is ruined by the 
existence of blow-up configurations. That kind of phe- 
nomenon is akin to the one found in simulations of lattice 
QCD, when the quark chemical potential is set around 
the threshold for fermion production After several 
atempts, we concluded it is unlikely that the numerical 
evolution can be stabilized by simple time step reduc- 
tions or the use of improved higher order Runge-Kutta 
schemes. However, from a purely mathematical perspec- 
tive, the source of instability in the Langevin simulation 
of the CDSG model is clear: it is related to the hyper- 
bolic functions defined in equations (|3.4I) . Their eventual 
fast growing, likely to happen close to phase transition 
points where field fluctuations are more intense, cannot 
be tamed by any finite precision numerical algorithm. 

Relying upon the idea of universality classes, we pro- 
pose a way out of the instability problem. One may argue 
that as far as the correlation length is much larger than 
the lattice spacing, the main interest is not the precise 
location, in the phase diagram, of the phases of the SG 
model. Instead, it is important to retain, in the contin- 
uum limit, the "topology" of the phase diagram and the 
correct singularity order of the transitions between dif- 
ferent phases. In loose words, we suggest to change the 
model, but not the large scale physics. We have found, 
in fact, that it is possible to deform the above set of 
equations, rendering them numerically stable, while pre- 
serving the BKT transition they are assumed to imply. 

To understand how the Langevin equations are to be 
modified, it is instructive to consider the following ordi- 
nary differential equation, directly inspired from the form 
of 



= — simp 



(3.5) 



Extending the field to the complex plane, by means of 
ip = ipi + iTp2, equation (|3.5|l becomes 



1p2 



sin -01 cosh -02 
cos ipi sinh -02 



The integral curves for this dynamical system, that is the 
vector field ("01, "02), are indicated in the top of Fig. 1. 
Due to the large "speeds" far from the tp2 = axis, nu- 
merical solutions exhibit bad uniform convergence prop- 
erties: for any choice of time step, a point close enough 
to one of the unstable manifolds would evolve towards a 
region of difficult numerical control. 




-202 
Re(v) 



3; 
E 



"MU'^" vH< Mi^ VM/ Mjr^ 



-2 

Re(v) 



(3.6) 



FIG. 1: The vector fields for the dimensionally-reduced dy- 
namical systems l|3.6^ and l|3.8^ . 

A simple modification of (|3.6|l leads to a dynamical sys- 
tem with uniform convergence properties, and the same 
asymptotic solutions. To get the alternative set of differ- 
ential equations, we perform, in l|3.5|l . the substitution: 

siiiijj = sin('0i -|- iV2) 
= cosh -02 sin ipi + i sinh 02 cos -01 
= cosh -02 (sin + i tanh ^02 cos V'l ) 
— > sin 01 + « tanh 02 cos "01 . 

Now, instead of H3.6|l . we will have 

ipi = -sin 01 , 

■02 = — cos V'l tanh -02 . (3. 



(3.7) 
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The integral curves for the above model are well-behaved, 
as it may be seen from the picture shown at the bot- 
tom of Fig. 1. Indeed, the vector field is bounded, for 
\\{ipi,ip2)\\^ = sin^(V'i) -l-cos^(V'i)tanh^(-02) < 1 and nu- 
merical convergence is never a problem in this framework. 
It is worth noting that the truncation prescription (|3.7(l 
keeps the original structure of the stable and unstable 
manifolds, as well as the positions of the fixed points 
in the (V'i,'!/'2) plane. Its role is basically to soften the 
vector field for the purpose of having stable numerical 
simulations. 

The same stabilization trick may be applied without 
any further complication to the CDSG Langevin equa- 
tions We get 



4^^ 



4^ 



4^ 



29? 0W _ 29i9o0k'^ 

5tanh(/30(2))cos(/30(i)) 

292^^1) + 2dido^ 

gsin(/30(i)) + 77L 

29i2<^i'^ + 29i9o<; 
gtanh(/30(2))(,Qg(^0(i)) ^ (3 9) 



4^ 



Numerical simulations were successfuly carried out for 
the above set of Langevin equations, using an elementary 
Euler evolution algorithm, for lattices with sizes 5x5, 
10x10, 15 X 15, and 20x20. Periodic boundary conditions 
were imposed on the field configurations. 
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FIG. 2: The mass order parameter $ is plotted as a function 
of l3 for different lattice sizes. The inset shows the data for 
the lattice with size 20 x 20. 

For each lattice size and a given value of /3, an ensemble 



of 6 simulations was considered, related to different seeds 
for the pseudorandom number generator. Each simula- 
tion consisted of 10^ iterations per lattice site. The SG 
coupling constant was set to g = 10.0. The time step was 
taken to be e = 2 x 10"'^. Averages were computed over 
a set of approximately decorrelated (in stochastic time) 
2000 configurations per simulation. 

In Fig. 2, the mass order parameter $ is plotted for 
(3 = 0.5, 1.0, . . . , 7.5. In more detail, we have evaluated, 
for a lattice of linear size L, 



$ = 



1 



2000 6 

6 X 2000 x 

k=l 1=1 



(3.10) 

where (fc, I) labels a given complex configuration of = 
4>{i,j) in the statistical sample. There is clearly a phase 
transition taking place around Pc = 2.5. The inset of 
Fig. 2 gives the best profile for the transition among our 
simulations, for the lattice with size 20 x 20. 




FIG. 3: The standard deviation analysis for the set of aver- 
aged evaluations of the mass order parameter <1>. 

In addition, we have carried out an estimate of the 
error bars related to the curves shown in Fig. 2. They 
are defined as 



(5$ 



\ 



where 



1 



2000 



2000 X L2 



k—l — l 



)l 



(3.11) 



(3.12) 



The results are shown in Fig. 3. As expected, fluctu- 
ations are larger around the phase transition, and get 
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smaller as the lattice size grows. For the lattice with size 
20 X 20, errors are just of the order of a few percent. 

Different lattice sizes were studied in order to estab- 
lish a finite size scaling analysis of the transition, re- 
ported in Fig. 4. Each one of the five straight lines 
draw in that picture corresponds to a fixed value of (3 
(/3 = 3.5, 4.5, 5.5, 6.5, 7.5), indicating that the scaling re- 
lation 



$ - L 



A 



(3.13) 



holds at the right side of the transition. According to the 
ideas of finite size scaling analysis |22], equation H3.13|l 
provides strong evidence for the fact that the whole re- 
gion (3 > 3.5 is a massless phase of the model under in- 
vestigation. No mass scale is involved; the only relevant 
scale turns to be the lattice linear size. Incidentally, as 
A ~ — 1 in the extended critical region, $ plays indeed 
the role of a moss order parameter, suggesting it is di- 
rectly related to the spectrum of massive excitations for 
f3 < (3c- We note that the phase transition separates a 
massive phase, with A = 0, from a line of critical points. 
This is the hallmark of the BKT phase transition, re- 
alized by the SG model. Therefore, we put forward the 
conjecture that the truncation prescription, implemented 
on the Langevin equations, is able to preserve the SG 
thermodynamical phases. 



Ml 
O 




Log(L) 



tion, revealed when one examines the results of simula- 
tions for the truncated CDSG model: while the scaling 
dimension of the mass order parameter depends quadrat- 
ically on f3 in the massless phase of the SG model, in the 
truncated context, it is approximately constant. That is 
not a surprise, however. Once we have deformed the orig- 
inal model, we cannot attribute the same scaling proper- 
ties to formally identical operators. Also, non-universal 
parameters, like (3c are not supposed to be the same in 
the alternative model. Actually, we have found /3c — 2.5, 
which is smaller than the value predicted for the SG 
model. 



IV. CONCLUSIONS 

The alternative, but equivalent, chirally decomposed 
version of the SG model yields an interesting stage for 
the test of numerical methods devoted to the analysis of 
systems governed by complex-valued actions. We have 
been able to show how a slight modification of the CDSG 
model at the level of the complex Langevin equations pro- 
vides a stable numerical framework, allowing us to iden- 
tify expected properties of the SG model, as the BKT 
transition. The essential idea of the truncation prescrip- 
tion, performed on the "bad behaved" hyperbolic terms 
contained in the complex Langevin equations, is to re- 
duce the degree of divergence of the unstable manifolds 
defined in the phase space, while preserving the thermo- 
dynamical phases of the original model. 

Since the CDSG model is similar to several rele- 
vant models found in condensed matter and high energy 
physics, it is likely that the rather elementary method 
of numerical stabilization addressed in this paper will be 
of utility in these contexts. Low-dimensional condensed 
matter bosonized models and the case of finite fermion 
density QCD are the ideal instances where an analogous 
line of numerical investigation could be applied. 

Furthermore, alongside with the numerical approach, 
a dynamical renormalization group analysis of the trun- 
cated stochastic differential equations is in order, to put 
the results established here on a firmer ground. In partic- 
ular, it is of fundamental importance to rederive the BKT 
renormalization group flow directly from the Langevin 
equations H3.9|l . a task which is under current investiga- 
tion. 



FIG. 4: Finite size scaling results obtained from the analysis 
of lattices with sizes 5 x 5, 10 x 10, 15 x 15 and 20 x 20. The 
inset shows the scaling dimension of the order parameter as 
a function of /3. 

An important difference to the SG physics calls atten- 
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